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Abstract 

Full three dimensional static and dynamic mean field calculations using collo- 
cation basis splines with a Skyrme type Hamiltonian are described. This program 
is developed to address the difficult theoretical challenges offered by exotic nu- 
clei. Ground state and deformation properties are calculated using static Hartree- 
Fock, Hartree-Fock+BCS and constrained Hartree-Fock models. Collective prop- 
erties, such as reaction rates and resonances, are described using a new alternate 
method for evaluating linear response theory, which is constructed directly on top 
of the static calculation. This provides a consistent description of the ground state, 
deformation and collective nuclear properties. Sample results are presented for the 
giant multiple resonances of 16 O. 



1 Introduction 



With the recent and future advent of new experimental facilities, such as 
the Radioactive Ion Beams facility at ORNL and GammaSphere, exploratory 
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detailed studies of many exotic nuclei will need correspondingly accurate the- 
oretical investigations. Examples of some of the properties of nuclei which 
will be studied include effective interactions, masses and ground state proper- 
ties, multi-dimensional energy surfaces and reaction rates. The results of these 
studies may also impact and address many astrophysical questions. 

The calculations discussed in this paper can be separated into two basic cat- 
egories. These are static mean field theories, such as Hartree-Fock, Hartree- 
Fock-Bogolyubov and constrained Hartree-Fock calculations; and using dy- 
namic linear response theory to study large-scale collective motion and reac- 
tion rates. 

Because of new capabilities in computer technology as well as innovations 
in numerical algorithms, a much greater level of sophistication can be incor- 
porated into numerical studies than in the past. In the studies presented in 
this paper, we incorporate the following ingredients. Collocation lattice repre- 
sentations using basis splines are used instead of traditional grid representa- 
tions or basis expansions. A full unrestricted three-dimensional calculation is 
performed, where no assumptions about spatial symmetries are needed. Spin 
degrees are included and one need not impose time reversal symmetry. As a 
result, to study nuclei with a large number of single particle states one requires 
grand challenge level computing resources. 

Collocation basis splines provide a superior representation than traditional 
grid techniques. The reason for this is that with basis splines, the gradient 
operator is represented as an operation upon a set of basis functions, which 
is much more accurate than the finite difference methods used in traditional 
calculations. At the same time collocation basis splines allow for the advanta- 
geous use of lattice representations. Here a full three-dimensional collocation 
lattice is used. 

It may be argued that there are existing accurate mean field calcualtions which 
use alternate basis expansions, such as sets of harmonic oscillator functions. 
One finds though that although these methods provide for accurate studies of 
more traditional nuclei, when one addresses exotic nuclei near the drip lines, 
such calculations begin to fail to accurately represent the physics. This can 
be traced to the fact that exotic nuclei, because of small binding energies, 
have very extended density distributions. One finds the harmonic oscillators 
bases typically used, due to their natural locality, are less able to represent 
accurately such extended objects than can a grid representation [|IJ. 

An alternate method has been developed to solve the linear response the- 
ory, which is a refinement of the technique used in ref. 0. Instead of solving 
a specific set of equations, the response of the system is derived using time- 
dependent Hartree-Fock [TDHF] theory. In this case a specific time-dependent 
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perturbing piece is added to the static hamiltonian to give H to t(t). The static 
Hartree-Fock solution is then time-evolved using H tot (t) in a TDHF calcula- 
tion. The Fourier transform of the result gives the response of the system at 
all energies. In this scheme one recovers both the response spectrum as well 
as the total transition probability amplitude corresponding to a given specific 
collective mode. 

The main advantage of this approach is that the dynamical linear response 
calculation is constructed directly on top of the static Hartree-Fock calculation 
and hence the static and dynamical calculations are calculated using the same 
hamiltonian description and the same degree of complexity. Therefore, there is 
a complete consistency between the static ground state of the system and the 
response calculations. One can then provide a coherent description of static 
properties of nuclei and of dynamical properties. This is important for example 
in f3— decay calculations of exotic nuclei, where reliable predictions are very 
sensitive to the deformation properties of the nucleus. Hence in this formalism 
consistent predictions of both the deformation and reaction rate properties 
are possible. 

Pairing correlations can be included in both the static and dynamic for- 
malisms. Presently, the BCS and Lipkin-Nogami prescriptions for including 
pairing have been used. Eventually we would like to use a self-consistent de- 
scription of pairing correlations, which would enable us to perform quasi- RPA 
studies of (3— decay for nuclei near the drip lines. With the recent advent of 
radioactive beam facilities, one would like to provide accurate and physical 
predictions of properties of exotic nuclei. 

Section 2 contains the theoretical discussion. Results are given in Section 3 
followed by a summary and conclusion in Section 4. 



2 Theoretical Discussion 



2.1 Continuous Equations 



The details of the derivation of the Hartree-Fock equations can be found in 
The result for a many-body Hamiltonian containing a one-body ki- 
netic energy and two- and three-body momentum dependent potential terms 
is a coupled set of non-linear partial differential eigenvalue equations, 



Xc 
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where Xa is a two- component vector (spinor) The effective nucleon-nucleon 
interaction is chosen to be of the zero-range Skyrme type, hence the Hamilto- 
nian h can be written in the following form (using natural units H — 1, c = 1, 
m = 1) 



h=~V 2 + W(p,rJJ), 

W = V N (r) + V c (r) , (2) 

where Vjq is the nuclear potential depending on various currents and densities 
and Vc is the Coulomb interaction. The densities and currents depend on the 
states Xa and are explicitly given by 



p(r) = J2^{\xUr)\ 2 + \x-(m (3) 

a 

r(r) = E^{ I (f ) | 2 + | Vx" (f ) | 2 } (4) 

a 

= Y J W a {Im[xi* {r)Vxt{r) + Xa*(r)VXa(r)} (5) 

a 

J{f) = -i J] «;«xS*(*0(V x a)x^'(f) . (6) 
Vc requires the solution of the Poisson equation in three-dimensional geometry, 

VV c (r) = -4vre 2 p(f) . (7) 



As can be seen from above, the solution of the system of equations (|1|) has 
to be obtained self-consistently and an accurate solution requires a good rep- 
resentation of various derivatives of the states \a- Currently, most HF and 
TDHF calculations are performed using finite difference lattice techniques. It 
is desirable to investigate higher-order interpolation methods which result in 
the improvement of the overall accuracy and reduction in the total number 
of lattice points. The lattice solution of differential equations on a discretized 
mesh of independent variables may be viewed to proceed in two steps: (1) 
Obtain a discrete representation of the functions and operators on the lattice. 
(2) Solve the resulting lattice equations using iterative techniques. Step (1) is 
an interpolation problem for which we could take advantage of the techniques 
developed using the spline functions The use of the spline collocation 

method leads to a matrix-vector representation on the collocation lattice with 
a metric describing the transformation properties of the collocation lattice. 
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2.2 Solution of the Discrete HF Equations 



The solution of the HF equations ([I]) is found by using the damped relaxation 



method described in Refs. |T0|JTT|| 



k+i 



o 



X k a -x B(E ) (h k -6 k a )x. 



where O stands for Gram-Schmidt orthonormalization. The preconditioning 



operator D is chosen to be [10,11 



D(£ ) 



E 



-If ml-l 

1+ E 



T 



where T denotes the kinetic energy operator. The solution is obtained through 
an iterative scheme as outlined below: 

(i) Guess a set of orthogonal single-particle states 

(ii) Compute the densities (6) 

(iii) Compute the Hartree-Fock potential 

(iv) Solve the Poisson equation 

(v) Perform a damping step (||) without orthogonalization 

(vi) Do a Gram-Schmidt orthogonalization of all states 

(vii) Repeat beginning at step (ii) until convergence 

In practical calculations we have used the damping scale value Xq ~ 0.05 and 
the energy cutoff E « 20.0. As a convergence criteria we have required the 
fluctuations in energy 



AE* = J(H*)-(H)'< 



(9) 



to be less than 10 -5 . This is a more stringent condition than the simple energy 
difference between two iterations, which is about 10~ 10 when the fluctuation 
accuracy is satisfied. The calculation of the HF Hamiltonian also requires the 
evaluation of the Coulomb contribution given by Eq. (^). Details of solving 
the Poisson equations using the splines are given in Ref. 



2.3 Collocation Basis Splines 



The static and time-dependent Hartree-Fock calculations are performed using 
a collocation spline basis in a three-dimensional lattice configuration. Basis 
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splines allow for the use of a lattice grid representation of the nucleus, which 
is much easier to use than alternate basis techniques, such as multi-dimensional 
harmonic oscillators. Also, for studies of exotic nuclei, because of weak binding, 
the density distributions tend to extend to large distances and hence one 
finds a large sensitivity to a harmonic oscillator basis, while for a lattice grid 
representation one needs to simply increase the size of the box. Traditional 
grid representations typically use finite difference techniques to represent the 
gradient operator. Collocation basis splines allow for the gradient operator 
to be represented by its action upon a basis function in a matrix form. This 
collocation method gives a much more accurate representation of the gradient, 
while maintaining the convenience of a lattice grid and hence provides a much 
more accurate calculation in the end. 

An Mth order spline function denoted by Bf 1 is constructed from piecewise 
continuous polynomials up to order M — 1. The set of points or knots {x{\ 
consists of the points where the spline functions are joined continuously up 
to the (M — 2) derivative. The basis spline functions have minimal support 
in that the zth spline functions is nonzero only in the interval (x^Xj+m), 
where the spline function, Bf 1 , is labelled by the first knot. For the case of 
using static boundary conditions for the space containing the N + 1 knots 
in one dimension, there must be M nonzero spline functions in each interval, 
hence N + 2M — 1 total spline functions make up the full basis, where M — 1 
spline functions extend beyond each boundary. One can also impose periodic 
boundary conditions, where there will then be N total spline functions. 

A function, f(x) continuous in the interval (x m i n , x max ) is expanded in terms 
of the spline basis functions: 

N+2M-1 

/(*)= E Bf{x)c\ (10) 



The expansion coefficients, d are derived from f(x) by evaluating f(x) at a 
specific set of points called collocation points, {x' a }. There are various ways 
of choosing the {x' a }. For odd order splines we have chosen the collocation 
points to lie at the center of each knot interval within the range {x m i n ,x max ). 



x a = , i = a, Va = 1, TV 



By evaluating eq. ( |10"D at {x' a }, a set of linear equations are constructed which 
constrain the coefficients, d: 

N+2M-1 

fWa)= E Bf\x' a ) C \ (11) 
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Since there are N+2M— 1 unknown coefficients in the case of static boundary 
conditions and N points x' a , it is necessary to introduce 2M — 1 additional con- 
straining equations as boundary conditions. Combining the functions Bf I {x' a ) 
and the boundary conditions into a square invertible matrix, B, the coeffi- 
cients, c % can be expressed as: 



w — » r 1 1 id 



'121 



where f a = f(x' a ) is the collocation representation of the function, f(x). 
Consider the action of an operator upon a function: 

N+2M-1 

Of(x)= £ [OBt\x)]c\ (13) 



If we evaluate the above expression at the collocation points and substitute in 
eq. fll2f ) for the d the following is obtained. 

N+2M-X ■„ 

of(x' a ) = df a = E K«)]E[b1 = Eq2/* ( 14 ) 

i 13 (3 



where now the quantity, Og is the collocation representation of the operator, 
O on the lattice. 

N+2M-1 .„ 

Oi= E [OB? (x' a )} [B-f (15) 



Note that the representation, Of, is not a sparse matrix. 

The function, f(x), and the operator, O, can both be represented on a lat- 
tice, i.e. the collocation points, through the use of the collocation basis spline 
method. This holds true for gradient operators, where the gradient of the ba- 
sis spline functions is required. The basis spline functions, B^(x) and their 

derivatives, 9 can be evaluated at the collocation points using iterative 

techniques. Through similar methods one can obtain the appropriate integra- 
tion weights: 

b 

1= f f(x)dx = J2w a f(x a ) } (16) 
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where the weights are given as : 



B™(x)dx 



B 1 



(17) 



For more details on the collocation basis spline method please see Ref. 



2.4 Linear Response Theory 



The linear response equations can be derived from a specific functional time- 
dependent perturbation of the TDHF equations |12| . To begin the proof a 
solution to the static Schrodinger equation is written as follows: 



H\^ s (0)) = E\M0))- (18) 

A time-dependent perturbing function is added to the static hamiltonian: 
H tot = H + H ex (t). (19) 

The external piece is defined as: 



H ex {t) = Ff{t) 

d 3 xn(x, t) F{x) 



fit), 



(20) 



where n(x, t) is the number density operator. F(x) corresponds to some par- 
ticular collective mode. The functions, F(x) and f(t) will be chosen later. 

At some time t = to the external piece of the Hamiltonian is turned on where 
\4>s(t)) is the solution to the time-dependent Schrodinger equation: 



= [H + H ex (t)] 



(21) 



Here the subscript V refers to the Schrodinger picture. The solution for the 
state vector can be written in the following iterative fashion. 

t 

\Ut)) = e-^ /7i |^(0)) - y- lSt/n J dt'HUt'MM) + - , (22) 

to 
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where the superscript '/' refers to the interaction picture. The expectation 
value of any operator, 0(t), is equal to 



(i> s (t)\Os(t)Mt)) = (&(0)|5 5 (0)|&(0)) + (23) 

t 

+m^jdt' [HUVAit)] |&(0)) + ... . 

to 

The linear approximation is made such that terms beyond first order in H[ x 
are neglected. 

If we define the operator 0(t) to be the number density operator, then using 
eq. (20) the fluctuation in the density becomes: 

5(n(x,t)) = ($M6s(t)\$.(t)) - (&(0)|d 5 (0)|&(0)> 

t 

= (h(0)\jrfMf d 3 x'F(x')f(t>) [njix'^M^t)] |&(0)> (24) 
to 

The retarded density correlation function is defined as: 

iD*(x, t; x', = B{t - ^M^™, (25 ) 

where hn = nn — (Rff) is the deviation of the number operator in the Heisen- 
berg picture. The density fluctuation can be written as: 

oo 

5(n(x,t)) = j- J dt' J d 3 x'D R (x,t;x',t')F(x')f(t'). (26) 



Using the Fourier representation of 9(t — t'), the Fourier transform of the 
density correlation function is: 



oo 

iD R (x, x'; uj)= J d(t - ty^-^iD^x, t; x', t') 



f {fpo\ns(x\ip n ){ipn\n s (x'\ipo) _ (^o\n s (x'\ij n ){ij n \n s (x \jj ) jj L 

where \ip n ) represents the full spectrum of excited many-body states of H . 
The Fourier transform of the density fluctuation then becomes: 
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5{n(x,uj)}= J dte tuJt 5(n(x,t)} 

— oo 

= d 3 x'D R (x,x';u)F(x')f(u), 

oo 

where f(u)= J dt'e lU}t ' f(t') 



(28) 



The linear response structure function, S(u>) is derived to be: 



f(u)S(u)= I d 3 x5(F\x)n(x,uj)) 

= ±-Jd 3 xJ d 3 x'F^(x)D R (x,x';uj)F(x')f(uj). 



(29) 



Combining eqs. (27) and (29) the imaginary part of the structure function 
then gives the total transition probability associated with F(x). 



IM[S(u)] = ~'£ 



d A x'{^ n \n x ^)\^)F{x' 



5 u- 



E n — E 



h 



E n > £ -(30) 



Note that this quantity is purely negative, where this feature can be used as 
a measure of the convergence of the solution. 

At this point instead of using the standard route of letting f(t) ^0 to recover 
the linear response equations, we choose an alternate technique to calculate the 
response. In this case we evolve the system in time and then fourier transform 
the result, where H ex (t) is a perturbing function. In this case f(t) is chosen 
to be a Gaussian of the following form. 



f(t)=ee-^)\ t >t 

f(u)=eJ—e-^ +vata , (31) 
V a 

where e is some small number (~ 10~ 5 ), chosen such that we are in the linear 
regime. The parameter, a, is set to be ~ 1.0 c/fm, which allows for a reasonable 
perturbation of collective energies up to ~ 150 MeV. The offset in time, t a is 
chosen such that the full gaussian is approximately within the time regime, 
t > t . 
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In practice to evaluate the TDHF equations the time-evolution operator is 
used to evolve the system. 



U(t,t )=T 



-i r df Htot (t') 



(32) 



where T [ ... } denotes time-ordering. Using infinitesimal time increments, the 
time-evolution operator is approximated by 



U(t 



n+l ) 



-| f t n+1 dt'H t „t(t') 
.±AtH tat (tn + ¥) 



fc! 



(33) 



N 

i + E 

fe=l 



/<)/, 



where the quantity B.\ ot is evaluated by repeated operations of H t 
wave functions. 



tot 



upon the 



The procedure is to then choose a particular form for F(x), using eq. (31) for 
f(t), and time evolve the system using eq. (33). The Fourier transform in time 
of the result then gives us f(u)S(u), from which the linear response structure 
function of the system is extracted. For a more detailed discussion please see 
Ref. pi. 



3 Results 



The static mean field calculation is analyzed in detail in Ref. [14|. For an ex- 



ample of some results using the static model in a HF+BCS calculation, please 



see Ref. [|L5] , where the results for various sulphur isotopes are presented. The 
pairing correlations are described by using a constant gap approximation of 
0.2 MeV. 

The static and dynamic Hartree-Fock calculations are performed using a 3- 
dimensional collocation lattice constructed with B-splines. Typically a 7th 
order B-spline is used in a (—10, +10 fm) 3 box with 20 3 — 24 3 grid points. 
The calculation can be performed with or without assuming time-reversal 
symmetry. The collective linear response may involve particle-hole interactions 
which are spin-dependent and not time-reversal symmetric. Therefore, for the 
correct collective content to be included one should not impose time-reversal 
symmetry in the linear response calculations |[16|| . In the results presented in 



this paper time-reversal symmetry is not imposed. In a comparison it is found 
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that imposing time-reversal symmetry causes small shifts in the position of 
the collective modes on the order of ~ 0.3 MeV. 

Calculations of isoscalar octupole; and isoscalar quadrupole collective modes 
are presented for 16 O using several parametrizations of the Skyrme interac- 
tion. Here the parametrizations known as skm* |T7J and sgll |Tj| are used for 



comparisons. A time step of At = 0.4 fm/c is used for the calculations. It 
is found that one can perform the time evolution for over 40000 time steps. 
The results shown here typically use 16384 time steps for a maximum time of 
« 6550 fm/c. 

Reasonable results are obtained if the parameter e in eq. (31) is chosen to fall 
in the range 2.0 x 10~ 4 < e < 2 x 10~ 7 . By varying the value of e, the amplitude 
of the time-dependent density fluctuation then scales proportionally to e, thus 
indicating that we are well within the linear regime of the theory. 

The linear response calculations require well converged initial static HF so- 
lutions. To test for the convergence of the static HF calculation the energy 
fluctuation, which is the variance of H as given in eq. @, is minimized. This 
measure of convergence tests directly how close the wave functions are to be- 
ing eigenstates of H and is independent of the iteration step size. For 16 it 
was found that static HF solutions with energy fluctuation less than about 
1.0 x 10~ 5 provided adequate starting points for the dynamic calculation, al- 
though the smaller the energy fluctuation the better. 

The dynamic calculations involve using eq. (33) to evolve the system. Since 
U(t, t') is an unitary operator, the ort honor mality of the system is preserved, 
therefore it is not necessary to re-orthogonalize the solutions after every time- 
step. The stability of the calculation is checked by testing the preservation of 
the norm of each wave function. The number of terms in the expansion of the 
exponent in eq. (33) is determined by requiring the norm to be preserved to 
a certain accuracy (typically to < 5.0 x 10~ 10 — > 1.0 x 10~ 8 ). 

The time-dependent perturbing part of the Hamiltonian is evaluated when the 
exponent in eq. (31) is greater than some small number, e cut . Since it is not 
difficult to evaluate the action of the external part of the Hamiltonian on the 
wave function, e cut is chosen to be very small, (1.0 x 10~ 10 ). To allow the fourier 
transform of f(t) to be evaluated easily, it is necessary to integrate t from — oo 
to oo and hence we would like the entire Gaussian of the perturbing function, 
f(t), to be included into the time evolution to the desired accuracy. The pa- 
rameter t a is therefore chosen such that the complete nonzero contribution of 



the time-dependent perturbation is included. t a = At ^2 + ^ 2 lo ^ ut ^ 



Pairing can be easily included using the BCS |T9| or Lipkin-Nogami |20 



con- 



stant pairing strength prescriptions. These two methods have been included 
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into the static Hartree-Fock calculations and can be easily incorporated into 
the dynamical calculation. For studies of ft— decay it will be necessary to 
include pairing, thus producing calculations of responses to quasi-RPA exci- 
tation modes. 

The computations were performed on massively parallel INTEL iPSC/860 and 
PARAGON supercomputers. 

For the study of the isoscalar quadrupole moment, the perturbing function 
F(x), introduced in eq. (20), is chosen to be the mass quadrupole moment, 
Q20 = 2z 2 — (x 2 + y 2 ). It turns out that other even multipole modes are also 
excited at the same time (i.e. Q i0 , Q 60 , ...). One can therefore study the effect 
of the coupling between the different excitation modes. The same holds true 
for the odd multipoles. 

In the top panel of Fig. 1 the time-dependent evaluation of the axial quadrupole 
moment, (Q2o{t)) is given for 16 using the skm* parametrization. 

(&o(t)> = (F(x)) = J d 3 x5(n(x,t))F^(x), (34) 

In this case the box is set to be (—10, +10 fm) 3 in size using 20 3 grid points. 
The periodic nature of the calculation is clear in this figure, where the smallest 
oscillation has a period of about 65 fm/c. 

Both a discrete Fourier transform [FT] and a continuous FT via Filon's method 
are used to calculate the FT of (Q2o(t)) to give (Q2o( u} )) — f( UJ )S2o(u). The 
frequency grid, Au is given by 2n/T max , so for a finer mesh one need only 
to time evolve to a larger maximum time, T max . The continuous FT is also 
limited by this criteria. The discrete and continuous FT's give essentially iden- 
tical spectrums, although the continuous FT recovers f{ui) more precisely and 
appears to give a more reasonable energy- weighted sum rule [EWSR] result. 
One can use the analytic expression for f(u) given in eq. (31) or a numerical 
result, where the difference in the spectrum is negligible. 

The lower panel shows the resulting response calculation. This quantity was 
derived in eq. (|3~0"D to be purely negative. The skm* results reflect this property 
quite well, where one sharp peak is seen at about 19.8 MeV. The calculation 
recovers about 92% of the exact EWSR, indicating excellent convergence. 

To test the sensitivity of the dynamic calculation to the box size, the calcula- 
tion shown in Fig. 1 is repeated using two larger boxes. The result in Fig. 1 
is compared in Fig. 2 to calculations using a box = (—11, +11 fm) 3 with 22 3 
grid points in the middle panel and one using a box = (—12, +12 fm) 3 with 
24 3 grid points in the lower panel. In all three cases the grid spacing is set 
to be 1 fm. One note is that since the continuum is represented by discrete 
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states in a box, the resonances are very sharply peaked. To best represent 
these peaks, T max is adjusted and chosen to give frequency grid locations as 
close to the peaks as possible. In the top panel it ended up that N = 16384, 
where T max = N x At, is an ideal choice. For the other two panels we found 
that choosing larger T max 's and hence larger iV's gave much better spectrums. 
One finds though that the spectrum changes as the box size is increased and 
in the middle panel there are 3 large peaks rather than just one. This disturb- 
ing result indicates that the spectrum as well as the peak positions appear 
to be significantly affected by the box size. We also find an effect when static 
boundary conditions are imposed instead of periodic ones. 

In Fig. 3 a calculation similar to Fig. 2 is shown, where the sgll Skyrme force 
parametrization is used. In this case the box also has a large effect on the 
calculated spectrum. The peak positions shift as well as the number of peaks. 
Also shown in the middle panel is a calculation which uses 22 3 grid points 
but with a denser grid. In this case there is also a significant change in the 
spectrum. 

The dependence observed in Figs. 2 and 3 can be traced to the representation 
of the continuum as given by the box. The lowest energy continuum states in 
the box correspond to the states with the fewest nodes. These states are of 
course dependent upon the size of the box. By increasing the size of the box, 
the lowest energy of these continuum states is reduced and the continuum is 
better represented in the calculation. Efforts are underway to increase the size 
of the box, by introducing unequal lattice spacing into the calculation. 

In Fig. 4 the isoscalar octupole response is shown for the skm* and sgll effective 
interactions. There are 3 sharp resonance structures at about 6 — 7 MeV, 
13 — 14 MeV and 18 MeV, where the sgll resonances are at the slightly higher 
energies. These three peaks are also found in a spherical RPA calculation 



using the same effective interactions ^TJ. The spherical calculation finds an 
additional broad peak for skm* and sgll centered at an energy of about 27 — 
28 MeV. This peak is weaker than the 3 peaks at lower energies and is not 
observed in the three-dimensional linear response calculation. In Fig. 4 one 
can see that the response is approximately purely negative, but that there are 
some violations of this feature. 



4 Summary and Conclusions 



A method for evaluating the linear response theory using TDHF is formally 
developed and implemented. This method allows one to construct the dynamic 
calculation directly on top of the static Hartree-Fock calculation. Therefore, by 
performing a sophisticated and accurate three dimensional static Hartree-Fock 
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calculation, we also have a corresponding consistent dynamic calculation. A 
coherent description of static ground state properties, such as binding energies 
and deformations, is given along with a description of the collective modes 
of nuclei. This feature is important in studies of (3— decay, where calculated 
reaction rates are very sensitive to the deformation properties of the nucleus. 

A collocation basis spline lattice representation is used, which allows for a 
much more accurate representation of the gradient operator and hence a cor- 
respondingly accurate overall calculation. Pairing correlations can be easily 
included into both the static and dynamic calculation, thus enabling us to 
look at quasi-RPA collective effects. 

Example calculations of 16 O are presented for the response functions corre- 
sponding to various isoscalar and isovector multipole moments. It is found 
that the time-evolution of the system in this full three dimensional calculation 
contains some dependences upon the box size, which need to be addressed. 
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Figure Captions 

Fig. 1 The fluctuations in the isoscalar axial quadrupole moment as a function of 
time are shown in the upper panel for 16 O using the skm* effective interac- 
tion for 16384 time steps. The size of the time step is 0.4 fm/c with a spatial 
grid of 20 3 in a cartesian box dimensioned (—10,-1-10 fm) 3 . The response 
for this calculation is shown in the lower panel. S(u) is obtained by fourier 
transforming the upper panel result and then dividing by f(u). 

Fig. 2 The lower panel of Fig. 1 is repeated here in the upper panel. The middle 
and lower panels correspond to the same calculation as shown Fig. 1 except 
the box sizes are (—11, +11 fm) 3 and (—12, +12 fm) 3 , respectively and the 
grid lattices are 22 3 and 24 3 , respectively. The number of time steps used, 
N are labeled on the graphs. 

Fig. 3 The same as Fig. 2, except the sgll parametrization of the Skyrme force is 
used and the number of time steps in the time evolution is different. Also 
shown in the middle panel as the dashed curve is a calculation using a 22 3 
lattice with a box size of (—10.6, +10.6 fm) 3 . 

Fig. 4 The imaginary part of the response function corresponding to the isoscalar 
octupole moment is shown for 16 and executed for 16384 time steps. The 
skm* and sgll parametrizations of the Skyrme force are used in the upper 
and lower panels, repectively with varying lattice grid and box sizes as 
labelled in the figure. 
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This figure "figl-l.png" is available in "png" format from: 



http://arXiv.org/ps/nucl-th/9410034vl 
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